This document reproduces the different steps followed before generating the shiny app.
library(tidyverse) # includes ggplot2 and readr
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag(): dplyr, stats
library(plotly) # interactive plots
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
Mset <- read_delim(file="http://www3.mbari.org/science/upper-ocean-systems/biological-oceanography/GlobalModes/Mset.txt",comment="%",delim="\t")
## Parsed with column specification:
## cols(
## YEAR = col_integer(),
## MONTH = col_character(),
## M1 = col_double(),
## M2 = col_double(),
## M3 = col_double(),
## M4 = col_double(),
## M5 = col_double(),
## M6 = col_double()
## )
Mset$time <- as.Date(paste(Mset$YEAR,Mset$MONTH,rep(15,length(Mset$YEAR))),format="%Y %m %d")
indices <- read_csv(file="http://www3.mbari.org/science/upper-ocean-systems/biological-oceanography/GlobalModes/climate_indices.csv",comment="%")
## Parsed with column specification:
## cols(
## year = col_integer(),
## month = col_integer(),
## MEI = col_double(),
## PDO = col_double(),
## NPGO = col_double(),
## EMI = col_double(),
## Nino3 = col_double(),
## Nino4 = col_double(),
## AMO = col_double(),
## NAO = col_double(),
## ATL3 = col_double()
## )
nb_time <- length(indices$year)
indices$time <- as.Date(paste(indices$year,indices$month,rep(15,nb_time)),format="%Y %m %d")
plot(indices$time,indices$MEI,type="l",xlab="Time",ylab="Index")
points(Mset$time,Mset$M1,type="l",col="red")
indexname1 <- "PDO"
indexname2 <- "NPGO"
iok1 <- !is.na(indices[,indexname1])
ipos1 <- indices[,indexname1]>0 & !is.na(indices[,indexname1])
ineg1 <- indices[,indexname1]<0 & !is.na(indices[,indexname1])
x <- indices$time
y0 <- rep(0,nb_time)
plot(x,t(indices[,indexname1]),type="l",xlab="Time",ylab=indexname1)
polygon(c(x[iok1],rev(x[iok1])),c(y0[iok1],rev(t(indices[iok1,indexname1]))),col="skyblue")
points(x,t(indices[,indexname2]),type="l",col="red")
Based on https://www.r-bloggers.com/r-single-plot-with-two-different-y-axes/
indexname1 <- "M1"
indexname2 <- "MEI"
iok1 <- !is.na(Mset[,indexname1])
x <- Mset$time
y0 <- rep(0,length(Mset$YEAR))
par(mar = c(5,5,2,5))
plot(x,t(Mset[,indexname1]),type="l",xlab="Time",ylab=indexname1)
par(new=T)
with(indices,plot(time,MEI,type="l",axes=F,xlab=NA,ylab=NA,col="red",pch=16,cex=1.2))
#plot(indices$time,t(indices[,indexname2]),type="l",axes=F,xlab=NA, ylab=NA,col="red")
axis(side = 4, col="red")
mtext(side = 4, line = 3, indexname2,col="red")
legend("topleft",legend=c(indexname1,indexname2),lty=c(1,1),col=c("black", "red"))
indexname1 <- "M3"
indexname2 <- "PDO"
max1 <- max(Mset[,indexname1],na.rm=TRUE)
max2 <- max(indices[,indexname2],na.rm=TRUE)
Mset[,"varpos"] <- Mset[,indexname1]
Mset[,"varneg"] <- Mset[,indexname1]
Mset$varpos[which(Mset$varpos > 0)] <- rep(0, length(which(Mset$varpos > 0)))
Mset$varneg[which(Mset$varneg < 0)] <- rep(0, length(which(Mset$varneg < 0)))
commonxM <- Mset$time %in% indices$time
commonxI <- indices$time %in% Mset$time
subsetM <- Mset[commonxM,indexname1]
subsetI <- indices[commonxI,indexname2]
iok1 <- !is.na(subsetM) & !is.na(subsetI)
r=cor(subsetM[iok1,],subsetI[iok1,])
legname <- paste(indexname2," (r = ",toString(floor(r*100)/100),")",sep="")
p <- ggplot() +
geom_area(data = Mset, aes(x=time, y=varpos), fill="blue") +
geom_area(data = Mset, aes(x=time, y=varneg), fill="red") +
geom_line(data = indices, aes(x=time, y=eval(as.name(indexname2))/(max2/max1), color=legname),size=0.5) +
scale_y_continuous(sec.axis = sec_axis(trans = ~.*(max2/max1),name = indexname2)) +
scale_x_date(limits = c(min(Mset$time),max(Mset$time)),expand=c(0,0)) +
labs(x ="Time", y = indexname1) +
theme(legend.position=c(0.15, 0.95),legend.key.width=unit(3,"line"),legend.background = element_rect(fill="transparent"),
legend.text=element_text(size=11)) +
scale_colour_manual("",
breaks = c(legname),
values = c("black"))
p2=plotly_build(p)
style(p2) %>% layout( legend = list(x = 0.01, y = 1) ) # needed because ggplotly removes the legend position